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Growing experimental evidence indicates that topological defects could serve as organizing centers 
in the morphogenesis of tissues. Here, we provide a quantitative explanation for this phenomenon, 


rooted in the buckling theory of deformable active polar liquid crystals. 


Using a combination 


of linear stability analysis and computational fluid dynamics, we demonstrate that active layers, 
such as confined cell monolayers, are unstable to the formation of protrusions in the presence of 
disclinations. The instability originates from an interplay between the focusing of the elastic forces, 
mediated by defects, and the renormalization of the system’s surface tension by the active flow. 
The posttransitional regime is also characterized by several complex morphodynamical processes, 
such as oscillatory deformations, droplet nucleation, and active turbulence. Our findings offer an 
explanation of recent observations on tissue morphogenesis and shed light on the dynamics of active 


surfaces in general. 


I. INTRODUCTION 


The development of multicellular organisms crucially 
hinges on internal and external mechanical cues which are 
transduced by the mechanosensing machinery of cells to 
give rise to system-wide spatiotemporal rearrangements 
and, eventually, to the formation of early embryonic fea- 
tures [1, 2]. These processes comprise a vast spectrum 
of active and passive forces, whose regulation relies not 
only on the molecular repertoire of tissue-forming cells 
but also on their shape, motility, and local organization. 
In the past decade, the conceptual framework of active 
matter [3] — namely, materials whose building blocks 
can autonomously move and perform mechanical work 
— has enormously contributed to classify the arsenal of 
physical mechanisms behind force generation and collec- 
tive migration in both eukaryotes and prokaryotes. How- 
ever, how these physical mechanisms can be harnessed 
to achieve biological functionality remains a major open 
question at the crossroad between developmental biology 
and mechanics [4-7]. 

One of the most fascinating hypotheses in this respect 
revolves around the possibility that multicellular systems 
could take advantage of topological mechanisms to cre- 
ate the conditions, in terms of reproducibility and ro- 
bustness, for the origin and maintenance of life. In par- 
ticular, defects have recently been identified by several 
in vitro studies as potential candidates for the role of 
“topological morphogens” in various biomimetic systems 
and model organisms. Keber et al. [8], for instance, 
demonstrated that protocells, consisting of a single layer 
of microtubule-kinesin active nematic liquid crystal en- 
closed in a lipid vesicle and powered by adenosine triphos- 
phate, relieve the mechanical stress, originating from 
the presence of four topologically required +1/2 discli- 
nations, by growing persistent tubular protrusions (Fig. 
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1A). More recently, the same mechanism has been in- 
voked in experimental [9-12] and theoretical [4, 13, 14] 
studies to explain the growth and regeneration of tenta- 
cles in Hydra (Fig. 1B). The accumulation of extensile 
stresses in proximity of integer and semi-integer discli- 
nations has also been proposed by Saw et al. [15] as a 
strategy to achieve cell apoptosis and extrusion in epithe- 
lial layers (see also Refs. [16, 17] for recent supporting 
evidence from simulations) and, more recently, by Guilla- 
mat et al. [18] as a route to the formation of multicellular 
protrusions in confined myoblasts (Fig. 1C). Long before 
the importance of defects for the development of non- 
planar features in tissues had been recognized, blister- 
like hemicysts — also known as “domes” — were already 
routinely observed in many epithelial cell cultures de- 
rived from renal cell lines (see e.g. Refs. [6, 19, 20]). 
An example of a dome, obtained from a monolayer of 
epithelial Madin-Darby Canine Kidney (MDCK) cells, 
was reproduced here and is shown in Fig. 1D (see also 
Movie S1 and S2 and Fig. S1). The central region of 
this dome, where the curvature is larger than elsewhere, 
features a multitude of topological defects (i.e. five- and 
seven-sided cells) with net positive topological charge: 
ie. icdome(6 — ci) > 0, with c; the number of sides 
of the ith cell. That is, the system exhibits an excess 
of positive disclinations (i.e. five-sided cells) where its 
curvature is maximal. This observation (see Secs. S1 
and S2 for further experimental evidence) suggests the 
possibility that, similar to the other examples illustrated 
above, the formation of cellular domes could be facili- 
tated by topological defects in combination with other 
system-specific mechanism, such as the injection of fluid 
under the cell layer which results in a focal detachment 
[21, 22]. 


Whereas these experimental studies have now convinc- 
ingly shown the existence of a correlation between topo- 
logical defects and certain morphogenic events, a clear 
theoretical picture of the mechanical nature of this be- 
havior is still missing. Here, we address this lack. Us- 


ing a combination of classical differential geometry, linear 
stability analysis and computational fluid dynamics, we 
demonstrate that confined active layers, here modeled as 
active polar liquid crystals, are unstable to the formation 
of protrusions in the presence of a +1 disclination. In 
extensile systems with positive flow alignment, the insta- 
bility originates from an interplay between the focusing 
of elastic forces, mediated by the defect, and a renormal- 
ization of the system’s surface tension by the active flow, 
arising in response to the perpetual injection of active 
stress under confinement. In the posttransitional phase, 
such competition leads to additional dynamical regimes, 
which includes oscillatory deformations, the nucleation 
of droplets, and a turbulent state with proliferating pro- 
trusions. By contrast, in contractile systems, the same 
phenomenon stabilizes the flat configuration, thus pre- 
venting the formation of protrusions. 

A precise account of all the examples of defect- 
mediated morphogenesis illustrated in the introduction 
requires, in principle, a great deal of biophysical details 
and a case-by-case approach. However, their very diver- 
sity suggests the existence of a general underlying mech- 
anism. Our goal is to identify and rationalize this mech- 
anism in the simplest possible setting, thus making it 
amenable for an in-depth analytical description. 


II. RESULTS 
II.1. The model 


Our model layer consists of a thin film of active polar or 
nematic liquid crystal, whose mid-surface M is spanned 
by the unit normal n and the pair of orthogonal tangent 
vectors gj, With i = 1, 2, and gij = gi - gj the associated 
metric tensor. The mean and Gaussian curvatures are 
denoted by H and K respectively (see Fig. 2A and, e.g., 
Ref. [23]). The local average orientation of the cells is 
described by the unit vector p = p'g;, with |p|? = ppi = 
1. The system free energy is given by 


> (Vip (V'P)| . 
(1) 
The first three terms on the right-hand side of Eq. (1) 
account for the energetic cost of deformations of the mid- 
surface, where y is the surface tension, Kp is the bending 
rigidity and «g is the Gaussian splay modulus [24]. The 
last term, on the other hand, is the Frank free energy [25] 
expressing the compliance of the system to a local distor- 
tion of the cellular polarization, with V; denoting covari- 
ant differentiation and «p denoting the rotational stiff- 
ness in the one-elastic-constant approximation. 
The steady state that we aim at describing comprises 
a dense population of cells freely translating and rotat- 
ing along a stationary mid-surface M. The former re- 
quirements amount to the following set of hydrodynamic 
equations for p and the momentum density pv = pv'g; 
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Figure 1. Examples of morphogenic events in the pres- 
ence of topological defects. (A) Active nematic protocell 
consisting of a monolayer of microtubules and kinesin enclosed 
in a lipid vesicle. The four topologically required +1/2 de- 
fects deform the lipid envelope, thereby creating tubular pro- 
trusions (adapted from Ref. [8]). (B) Example of Hydra fea- 
turing +1 disclinations in proximity of the mouth, the foot 
and the tip of each tentacle, and two —1/2 defects at the base 
of each tentacle (adapted from Ref. [11]). (C) Multicellular 
protrusion in layers of collectively migrating myoblasts, fea- 
turing a +1 disclination at the tip (adapted from Ref. [18]). 
(D) Example of a dome obtained from an initially flat mono- 
layer of cultured MDCK GII cells. The central region of the 
dome, where the curvature is larger than elsewhere, hosts sev- 
eral topological defects (i.e. five and seven-sided cells), with 
a net positive topological charge. On the left-hand side a 
cross section and a top view of the dome with a superim- 
posed tessellation are shown (red, F-actin; green, E-cadherin; 
blue, nuclei). 


on the tangent plane [26-29]: 


P(e + oY) = Vy (of? +08) + f, (2a) 


(3 + vV k) = (9 — p'p’) (Aujnp” = wep” + Ph;) ; 
(2b) 


The total cellular mass M = fdAp is assumed to be 
constant in time, although the number of cells could in 
principle fluctuate as a consequence of cell division and 
apoptosis. Thus p = const and V,;v’ = 0. This as- 
sumption is justified by the fact that, in cellular sys- 
tems, the speed v of migrating cells is orders of mag- 
nitude lower than the speed of sound cs. Because the 
typical magnitude of density fluctuations is 6p ~ Ma’, 
with Ma = v/c, the Mach number, this implies, to first 
approximation, dp ~ 0. The three terms on the right- 
hand side of Eq. (2a) correspond, respectively, to the 


hydrodynamic stress tensor o? = —P g” + 2nu’ 
where Ph is the pressure, 7 is the shear viscosity and 
uti = (Vivi + Vv") /2 is the strain-rate tensor— the ac- 
tive stress ot and the force per unit length fi, originating 
from the distortion induced by defects and Gaussian cur- 
vature and driving the so-called elastic “backflow” [30]. 
The latter can be derived from expressing the polariza- 
tion gradients in terms of the geometric potential x, sub- 
ject to the Poisson equation V?x = K — pg, where pa 
is the topological charge density [29, 31-33]. Explicitly, 
fi = —2KrpaV'x (see SI). Thus, at equilibrium, pa ~ K 
so that defects place themselves in regions of like-sign 
Gaussian curvature to minimize the elastic free energy 
[33]. Last, the active stress tensor embodies the contri- 
bution of the force dipoles autonomously generated by 
the cells and, following a standard construction, can be 
expressed as otf = a(p'p! — g” /2) [3, 34-37]. The phe- 
nomenological constant a quantifies the magnitude of the 
cellular forces and is positive (negative) for contractile 
(extensile) systems. 

Equgation (2b) governs the rotational dynamics of the 
cells which, in turn, is dictated by the coupling with 
the local flow field, expressed by the first two terms 
on the right-hand side of the equation, with » the flow 


alignment parameter and w = (Vivi — VJv')/2 the 
vorticity tensor. The last term is the molecular field 
hi = —bF/ép’ = krV*p; which drives the system to- 


ward the ground state of the free energy, with T7! the 
rotational viscosity (see e.g. Ref. [25]). 

Last, the stationarity of the mid-surface M requires 
the net force acting along the normal direction n to van- 
ish, hence 


Ki; (oi + ci!) + fet feso, (3) 


where K;; is the curvature tensor (see Fig. 2A). The nor- 
mal forces per unit length fẸ and f?’ on the left-hand side 
of Eq. (3) are found from minimizing the free energy F 
(see SI) and are given by f? = 27H + 2kgH(K — H?) — 
kpV2H and f} = 2np(2Hg — KY)V;V jx + 2kp(KY — 
Hg” )VixV;x [29, 38-41]. We stress that Eq. (3) re- 
duces to the Helfrich shape equation, fẹ = 0, in the limit 
a = kp = 0 [23, 38], and to the van Kármán equation, 
fg + 27H = 0, in the limit a = kg = 0 [42, 43]. 


II.2. Stability of defective active layers: Topology, 
geometry and morphogenesis 


In this section, we investigate how topological defects 
and activity conspire to render an initially flat layer un- 
stable to the growth of protrusions. While Eqs. (2) and 
(3) hold for arbitrary conformations of the midsurface 
M, here we consider the simpler case of an axisymmetric 
surface and a liquid crystal with polar symmetry, where 
the local polarization and velocity fields depend solely 
on the distance r < R from the symmetry axis, with R 
the radius of the layer: i.e., p = p(r) and v = v(r). 


Under this assumption, the polarization field inevitably 
features a +1 topological defect at the center of the sys- 
tem. Thus, taking pq = 270(r), with 6(-) a delta function 
on M, and setting y = 0 at r = R without loss of gener- 
ality, the geometric potential can readily be found to be 
x = —log(r/R). 

To make progress, we parameterize the midsurface M 
in a Monge patch, so that the position r = r(r,y) of a 
generic point is given by r = rg, + he,, with e, a unit 
vector in the z direction and h = h(r) the height. Follow- 
ing a standard approach of membrane physics (see, e.g., 
Ref. [23]), we then assume that |Vh| < 1, so that one 
can ignore terms of order O (|Vh|?) and higher in Eqs. 
(2) and (3). Under such a small-gradient approximation, 
gij © Oj, H ~ -V7h/2, K ~ 0, and Eqs. (2) reduce 
to their flat counterpart. Moreover, axial symmetry and 
incompressibility yield v” = 0, so that v = 1/rv%gy, 
with v? = v*(r), and p = coseg, + 1/r sinegy, with 
€ a constant and (r,y) polar coordinates on the h = 0 
plane. Solving Eq. (2b), we find e = tarccos(—1/A)/2 
under the assumption that |A| > 1. Similarly, neglecting 
unimportant inertial terms [44] and solving Eq. (2a) for 
no-slip boundary conditions gives, after standard alge- 
braic manipulations (see e.g. Ref. [45-47]): 


ITZ 2 
v? = 72 ESA r log 2 7 (4a) 
n 2 R 

a 
A 


consistent with a classic result by Kruse et al. |45]. It is 
worth stressing that, by virtue of Eq. (4b), the pressure 
nominally diverges at the origin, where the +1 defect is 
located. Similar to other divergences stemming from the 
continuous description of defects, however, this one too 
can be regularized by introducing the coherence length 
le, setting the short length scale cut-off at which the 
continuous theory breaks down, because the magnitude 
of the polarization field, here assumed to be constant, 
drops in proximity of a defect. Thus e < r < R. Last, 
replacing Eqs. (4) in Eq. (3), gives a stability condition 
for the flat conformation of the midsurface M in the form 
of a fourth-order linear partial differential equation with 
position-dependent coefficients: 


E 
Pi = log p “FTPa, (4b) 


fB Yth — yet V?h + n o, (2>) =0, (5) 

2 r r 
where Yer = y — Ph and kef = Kp + ar?/(2A) are the 
effective surface tension and Frank’s elastic constant, re- 
spectively, which include the renormalization resulting 
from the active flow. Note that the velocity field does 
not enter Eq. (5) explicitly, as the term Kj;u’? vanishes 
identically. Moreover, M is clamped at the boundary 
such that h and all its derivatives vanish at r = R. We 
stress that the assumption of rotational symmetry re- 
stricts us to the low-activity regime. Upon increasing the 
activity, the system would eventually enter a turbulent 
regime, thereby losing its spatial symmetry. Evidently, 
this regime cannot be captured by our solution. 
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Figure 2. Post-buckling geometry of layer. (A) Sketch of 
the midsurface M embedded in R?. In the table, the cross and 
the dot refer to the vector and scalar product, respectively, 
with respect to the Euclidean metric of R®. In the figure and 
throughout the article bold letters indicate R? vectors and 
Latin indices indicate surface coordinates on M. The coor- 
dinate system on M is denoted by {g:,n} and it defines the 
metric gij and the curvature tensor Ki; of the surface, hence 
the mean curvature H and Gaussian curvature K. (B) In the 
top middle a sketch of the buckled surface with a +1 defect at 
the center. The disk on the top right is a view from the top. 
The black arrows are the azimuthal velocity field, the color 
scale is its magnitude, and the white arrows are the director 
field. In addition, we solved Eqs. (2) and (3) numerically 
to plot the height profile h(r/R) of the surface for different 
values of the active number A = al /«r. Comparing the ac- 
tive with the passive (A = 0) case, we see how an extensile 
activity (A < 0) favors the buckling, while a contractile one 
activity inhibits it. The values of the constants used are the 
following: y = 0.012, T = 1, A= 1.1, kr = 0.02, kg = 0.01, 
and 7 = 5/3. 


Now, to gain insight into the stability of the system 
with respect to the formation of protrusions, we sepa- 
rately consider three limit cases, namely, (i) kg = 0 and 
a = 0; (it) kg = 0 and a Æ 0, and (iii) kg # 0 and 
a #0. From a physical perspective, the first two cases 
correspond to layers of passive and active subunits, re- 
spectively, whose elastic stiffness is sufficiently small to 
be considered negligible, but where surface tension pe- 
nalizes an increase in the area. The vanishing bending 
modulus, in particular, renders highly curved features, 


such as kinks or cusps, energetically acceptable. The first 
scenario has been considered by Frank and Kardar [48] 
and corresponds to the case of a passive liquid-crystalline 
disk-like interface plagued by a +1 disclination in the 
center. In this case Eq. (5) can be integrated exactly, to 
give 


(r? — RÈ) O-h=0, (6) 


with RŽ? = kp/y. Because of our boundary conditions 
for h, this equation admits a nontrivial solution only if 
R > R. = Vkr/y. Even Eq. (3) is exactly integrable 
in this case and M is found to be a parabolic pseudo- 
sphere —namely a surface of constant negative Gaussian 
curvature— whose height and mean curvature diverge as 
r — 0 [48]. Similarly, in the second scenario, Eq. (5) 
reduces to Eq. (6), but with Re given by the solution 
of the transcendental equation Kp /R? = y — a/(2X) — 
a/Alog R/R.. As in the previous case, we find that there 
is a non-trivial solution only if R > Re, thus the flat con- 
formation becomes unstable to the growth of protrusions 
when the radius R of the active layer exceeds the critical 
radius 

R= |. (7) 

Y ox 


Equivalently, for fixed R values, Eq. (7) provides a 
threshold in |a|, above which the layer becomes unsta- 
ble to buckling. Note that for A > 1 the activity reduces 
(increases) the effective surface tension in the presence of 
contractile (extensile) stresses and vice versa for negative 
A (see also Ref. [41]). Conversely, if y < a/(2A), then 
Eq. (5) has no nontrivial solutions and the flat solution 
is always stable. 

Last, in the third and most generic scenario, Eq. (5) 
cannot be reduced to Eq. (6) and one cannot find a stabil- 
ity criterion as simple as that expressed by Eq. (7). How- 
ever, as for the passive instability discussed in Ref. [48], 
we expect the bending stiffness xg to merely renormalize 
Frank’s elastic constant kp. In addition, a finite bending 
rigidity makes cusps, kinks and other singularities charac- 
terized by a diverging curvature energetically prohibitive; 
hence we expect the tip of the protrusion to become in- 
creasingly smooth with increasing Kg. To substantiate 
the latter statement, we have lifted the small-gradient ap- 
proximation and numerically integrated Eqs. (2) and (3) 
on an axisymmetric surface. The equations do not decou- 
ple in this limit, therefore it is not possible to condense 
them in a single equation for the height h. However, 
one finds from the incompressibility equation and with 
p = 1/g coseg, + 1/r sineg,, where g = y1 + (0,-h)?, 
that v = 1/r v?gą, and the angle € as well as the pressure 
P, are left unchanged. The remaining fields, v? and h, 
can be numerically computed and h is shown in Fig. 2B 
for different values of activity. Note that the solutions 
found for v? are essentially identical to the analytical 
solutions obtained using the small-gradient expansion. 

In summary, the presence of a +1 defect renders the 
planar configuration of an active polar layer unstable to 
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Figure 3. Dynamical regimes in active layers. (A to E) Defect-mediated buckling and activity-driven shape deformations 
of an active polar layer obtained from numerical simulations (see Sec. II.3). The corresponding phase diagram is shown in (F). 
In all simulations we fix the dimensionless surface tension © = 0.6 and vary the dimensionless activity A. The white vectors 
mark the polarization field, while the color code in (G) refers to the local magnitude of the flow. At small activity (yellow 
pentagon in the phase diagram), we find a vortical flow structure and the interface is flat outside the defect core. This is 
illustrated in (A) where A = —0.2 and the inset shows the spiral pattern in the polarization field. Beyond the threshold of the 
buckling instability (green squares), the layer attains a funnel-like profile with the topological defect at the end of the hollow; 
this is shown in (B). At slightly higher activity, just below the transition to the region of droplet production (A = —0.4), 
we find that the interface begins to oscillates between this funnel-like state and the configuration shown in (C). Next, in (D) 
we show snapshots taken at different times during the process of droplet nucleation at A = —0.55 (cyan circles). Last, in 
(E), we show a chaotic configuration of the layer in the active turbulent regime (A = —1.0, blue diamonds). (H) The time 
evolution of the normalized distance in the L2-space between the simulated profile of the phase field ¢ and the equilibrium one, 
ġ* = tanh (z/€), with £ the equilibrium interface width between the two isotropic fluids and z the coordinate in the direction 
normal to the interface. The inset on the top left in (H) refers to the case at A = —0.4 and shows the periodic oscillations 
between the configurations of (B) and (C) for A = —0.4 in terms of the normalized L2-distance. 


buckling. The instability arises from the fact that the de- previous section, is a notable example of how the inter- 
fect introduces an angular deficit in the orientation of the play between the topology of the polarization field and 
polarization field that a like-sign Gaussian curvature can the geometry of the flow cooperatively render an ini- 
compensate, thereby mitigating the distortion sourced by tially flat layer unstable to the growth of protrusions. 
the defect [33]. Despite having a passive origin, this in- In this section, we look at the evolution of the instability 
stability is affected by the hydrodynamic flow fueled by to unveil how topological defects influence the morphol- 
the layer’s extensile or contractile activity. This flow ogy of the protrusion in the post-buckling scenario. To 
changes the later pressure acting on an arbitrary fluid this end, we numerically simulate an active polar layer, 
patch, resulting in an additional compression (for exten- whose thickness € is comparable in magnitude with the 
sile activity) or stretching (for contractile activity) of the coherence length ¢., which, in turn, represents the short- 
layer, which, in turn, inevitably interferes with how the est length scale in our continuum description. The ac- 
layer itself deforms out of plane, thus making the system tive layer is sandwiched between two Newtonian fluids, 
respectively more or less probe to buckling. whose relative concentration is described by the phase 
field 6 = g(r) [49, 50]. In the limit € — 0, the simulated 
interface can be factually interpreted as the mid-surface 
M of Sec. II.1 [18, 51]. A description of the model is 
provided in the Materials and Methods section and SI. 
We numerically integrate the three-dimensional (3D) hy- 
drodynamic equations of this diffuse interface model in 
a cylindrical container with homeotropic boundary con- 


II.3. Buckling, oscillatory regime and droplet 
nucleation 


The interplay between defect-mediated stress focusing 
and the elasticity of the active layer, illustrated in the 


ditions along the base, so that the equilibrium configu- 
ration in the absence of activity is stationary and char- 
acterized by a +1 disclination at the center of the disk. 
The interface is slightly deformed at the defect position 
to accommodate the tendency of the liquid crystal to es- 
cape in the normal direction. Outside the defect core 
Lle the interface is flat, in agreement with our analytical 
prediction in Sec. II.2, since, in the absence of activity, 
we have buckling only in a disk of radius Re = \/5/3 le- 
In the following we will restrict our analysis to extensile 
systems as we find that contractile ones do inhibit a buck- 
ling instability, in agreement with what is predicted by 
our analytical argument. The main control parameters 
that we varied in our numerical experiments are (i) the 
dimensionless activity A = af?/Kp and (ii) the reduced 
surface tension © = yé./Kp. 

As the dimensionless activity A is increased, four dif- 
ferent regimes are encountered. At small A values, the 
surface is flat, whereas the polarization evolves into a spi- 
ral configuration coupled to a vortex flow (Fig. 3A), in 
agreement with our prediction in Sec. II.2 and Ref. [45]. 
As activity is further increased, the interface undergoes 
a buckling instability (light green squares in Fig. 3F), 
which results in the development of a protrusion, fea- 
turing a +1 defect at the tip. Thanks to our computa- 
tional approach, we notice that once the instability is set 
by the previously described competition between defect- 
mediated stress focusing and elasticity, the development 
of nonplanar features is further facilitated by the length- 
ening of the vortical flow field out of plane, a process 
called vortex stretching, which occurs in 3D fluids as the 
result of conservation of angular momentum [52]. The 
transition line between the flat and buckled configura- 
tions in the phase diagram of Fig. 3F shows a linear be- 
havior in the A — © plane, consistent with the analytical 
criterion expressed by Eq. (7). 

The competition between active and elastic stresses 
may eventually give rise to large perturbations around 
the stationary cuspidal profile, with the surface oscillat- 
ing from the singular configuration with negative curva- 
ture of Fig. 3B to the smooth configuration of Fig. 3C. 
The dynamic of the oscillations can be captured by mea- 
suring the temporal evolution of the distance of the sim- 
ulated interface from the flat configuration, shown in 
Fig. 3H, and its inset for the case at A = —0.4 and 
£ = 0.6 (see also Movie $3). The frequency of the oscil- 
lation linearly increases with |a| along the transition line, 
which can be rationalized either by dimensional analysis 
or from the balance of non-equilibrium stresses and vis- 
cous ones. This oscillatory instability lies at the heart 
of another notable behavior which is observed when ac- 
tivity is further increased. In this case, the surface’s en- 
tropic elasticity is no longer able to counteract the active 
stresses fueling the growth of the protrusion (top panels 
in Fig. 3D and Movie S4). This results in a steady in- 
crease of the passive stresses along the protrusion, which 
eventually results in the breaking of the interface and 
the consequent nucleation of a droplet (bottom panels in 


Fig. 3D and Movie $4). The emulsified droplets are there- 
fore enclosed in a thin active polar shell, which separates 
the interior of the droplet from the outer Newtonian fluid. 
For each droplet, the polarization field develops two boo- 
jums (i.e. +1 defects) as required by the Gauss-Bonnet 
theorem. 

In the limit of very large activity (blue diamonds 
in Fig. 3F), the dynamics becomes fully chaotic, a 
regime which is known as active turbulence in the lit- 
erature [26, 44, 53, 54]. The active layer is characterized 
by the proliferation of many amorphous protrusions (see 
Fig. 3E and Movie S5) which elongate under the straining 
effect of bending deformations in the polarization pat- 
tern. A similar chaotic state is also observed for very 
large contractile activity. 

Last, we stress that the results in Sec. II.2 are quanti- 
tatively unchanged by the frictional interaction between 
active layer and the surrounding Newtonian fluid, which 
the diffused-interface model summarized in this Section 
naturally accounts for. This can be understood by the 
fact that, at low Reynolds number, such an interaction 
gives rise to a damping force fe = —Cv, with Ç a drag co- 
efficient, in the active layer [55]. This, in turn, introduces 
an additional length scale, i.e. l = /n/¢, reflecting 
the competition between viscous (i.e. momentum con- 
serving) and frictional dissipation. That is, for distances 
smaller than @¢, friction does not produce substantial ef- 
fects and momentum is approximately conserved as in 
the absence of friction. By contrast, at distances larger 
than ¢¢, frictional dissipations takes over and the momen- 
tum density decays exponentially. Because the defect- 
mediated buckling transition discussed here is highly lo- 
calized around the defect core, friction will not produce 
any substantial change, unless the coefficient is unrealis- 
tically large. 


III. DISCUSSION AND CONCLUSIONS 


Here, we elucidated how the interplay between topol- 
ogy, geometry and hydrodynamics enables the develop- 
ment of nonplanar features in active liquid crystals, as re- 
cently observed experiments on biological and biomimetic 
systems and, more prominently, in eukaryotic cell layers 
(Fig. 1 and Refs. [9-12]). Whereas the amount of bio- 
chemical detail necessary for a thorough account of all 
aspects of this phenomenon often results in a complex 
tangle that hinders fundamental understanding, here, we 
followed a more generic approach, by focusing on the dy- 
namics at the mesoscopic scale and retaining only the 
orientational and elastic degrees of freedom, which all re- 
alizations of defect-mediated morphogenesis have in com- 
mon. 

By looking at an axisymmetric surface plagued in the 
center by a +1 disclination, we analytically demonstrated 
that defective active layers with liquid crystalline order 
are unstable to out-of-plane deformations and obtained 
a stability criterion in terms of system size, orientational 


stiffness, surface tension and active stresses. This mech- 
anism allows defects to serve as topological morphogenes 
for the formation on nonplanar features, such as domes 
and protrusion. Such an instability originates from the 
competition between the focusing of the elastic forces, 
mediated by defects, and the renormalization of the sys- 
tem’s surface tension by the active flow. Upon modeling 
the cell layer as a diffuse active polar interface and turn- 
ing to computational fluid dynamics, we further inves- 
tigated the posttransitional regime and constructed an 
exhaustive phase diagram. At low activity, we recover 
the results of our analytical approach. Upon increasing 
the activity, we first find a regime where the layer os- 
cillates periodically between two different buckled states 
and then a regime where the high activity breaks up the 
interface leading to continuous droplet nucleation. Last, 
at very large activity, we see a turbulent regime which is 
characterized by the chaotic proliferation of protrusions. 


Whereas undoubtedly simplified with respect to the 
seaming endless complexity of living matter, our study 
contributes to shed light on how developing organisms 
can take advantage of topology to achieve biological 
organization from physical mechanisms, with potential 
application to various biomechanical processes such as 
morphogenesis, embryogenesis, cancerogenesis and vesi- 
cle formation. As mentioned in the introduction, +1 
disclinations are not the only type of defects encountered 
in morphogenesis, but nematic (i.e. +1/2; see Fig. 1A,B) 
and, in general, p—atic defects (e.g. +1/6, see Fig. 1D), 
are equally common and we will investigate them in the 
future. Last, while there are several experiments observ- 
ing the role defects play in the formation of protrusions, 
we are not aware of any study investigating the other 
regimes, particularly droplet nucleation. Naturally, the 
question occurs whether these other regimes also play a 
role in biological systems. 


MATERIALS AND METHODS 


Experimental setup 


Parental MDCK GII cells (provided by M. Gloerich, 
UMC Utrecht) were grown on non-coated coverslips un- 
til tissue buckling. F-actin, E-cadherin and nuclei were 
stained on fixed samples. Samples were imaged at high 
resolution on a spinning disk confocal microscopy setup. 
Cell boundaries were identified from a maximum inten- 
sity projection of a z-stack of the top part of the dome. 
Fiji software was used for the orthogonal view. 3D re- 
constructions were done by Imaris Viewer 9.7.0 and z di- 
rections were corrected for spherical aberration and axial 
distortion [56]. See the SI for details. 


Numerical simulations 


The dynamical fields of the model used for simulations 
are the incompressible flow field v = v(r,t) (V -v = 0) 
and the polarization field P = P(r,t), which is confined 
at the interface between two isotropic fluids, whose rel- 
ative concentration is encoded in the scalar phase field 
$ = (r,t). Note that in this section P and v are three- 
dimensional vectors, and V now denotes differentiation in 
R°. The equilibrium is defined by a generalized Landau- 
de Gennes functional [49, 50]: 


_ a 2 Q 4, ko 2 
F(6.P\= fav -56+ 30+ 2 (v0) 
+4 ($ IPP | 7lPI') T VPI? + Ë (P-V)? l 


(8) 


where the constants a and kg are model parameters re- 
lated to the surface tension and the width of the inter- 
face by y = \/8akg/9 and € = y/2k¢/a, respectively. 
To confine the polar liquid crystal at the ¢-interface, the 
parameter 7 = u(V@) is chosen to be a function of Vo 
such that 7 = —1 if |V¢@| is larger than a suitable thresh- 
old and 0 otherwise. In addition, the coupling between P 
and V¢ ensures tangential anchoring of the liquid crystal 
for 6 > 0, so that the polarization field actually lays on 
the interface. The bulk constant Ag fixes the coherence 
length of the liquid crystal le = \/«r/Ao, which con- 
trols how fast the order parameter P drops to zero from 
its equilibrium value |P| = 1 in proximity of a topologi- 
cal defect. Note that in the limit of vanishing interfacial 
width, € — 0, such a phase-field model can be mapped on 
the analytical model of Sec. II.1 (see e.g. Refs. [13, 51]). 
The dynamics of the system is governed by the following 
set of partial differential equations: 


(3 +v: Vv = V : (op +a) , (9a) 
(3 +v- V)P =w. P+\u.P+rh, (9b) 
AO+V- (dv) =V- (uv) ; (9c) 


Equation (9a) is the Navier-Stokes equation which rules 
the hydrodynamics of the system in the full three- 
dimensional space. Here, the stress tensor has been di- 
vided in a passive and an active part. The former in- 
cludes, in turn, three terms: a hydrodynamic, an elastic 
and a phase-field contribution op = Oh + Ce + oj, whose 
explicit expressions are given by 


Oh = —P,1 + n [Vv + (Vv)"] j 


A 
2 


) 1—kyV¢V¢—-BPV¢. 
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o 
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Here, P, is the hydrodynamic pressure, 7 is the shear 
viscosity, À is the flow-alignment parameter, and I~! 
is the rotational viscosity. f denotes the free energy 
density, such that F = [dV f, h = -F/P is the 
molecular field and u = [(Vv) + (Vv)™]/2 and w = 
[((Vv) —(Vv)*]/2 are the strain rate and vorticity tensor, 
respectively. The three-dimensional active stress tensor 
Oa is given by 


oO, =a (pp = : IPn) l (10) 


Last, Eq. (9b) is the Ericksen-Leslie equation for the 
polarization field P in three dimensions, whereas the 
dynamics of the concentration field is governed by the 
advection-diffusion equation, Eq. (9c), with u the mobil- 
ity parameter. 

The dynamical equations have been integrated by 
means of a hybrid lattice Boltzmann (LB) method [57], 
where the hydrodynamics is solved through a predictor- 
corrector LB algorithm, while the dynamics of the or- 
der parameter has been treated with a finite-difference 
approach implementing a first-order upwind scheme and 
fourth-order accurate stencil for the computation of spa- 
cial derivatives. The hydrodynamics was integrated in a 
cylindrical geometry with radius R = 324. under soft 
homeotropic boundary conditions for the polarization 
field, while the concentration field satisfies neutral wet- 
ting boundary conditions at the top and bottom walls, 
respectively, where 6 = +1. The system is initialized 
with ¢ = 1 (¢ = —1) in the top (bottom) half of the sys- 
tem and the polarization field in an aster configuration at 
the interface in the center of the system and 0 otherwise. 

The numerical code has been parallelized by means of 
Message Passage Interface (MPI) by dividing the compu- 
tational domain in slices and by implementing the ghost- 
cell method to compute derivatives on the boundary of 
the computational subdomains. Runs have been per- 
formed using 64 CPUs on a computational box size of 
96 by 96 by 256, for at least 10° LB iterations (corre- 
sponding to ~ 35 days of CPU time on Intel Xeon 8160 
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I. EXPERIMENTAL DETAILS. 


Cell culture. Parental Madin-Darby Canine Kidney (MDCK) GII cells (kindly provided by M. Gloerich, UMC 
Utrecht) were cultured in a 1:1 ratio of low glucose DMEM (D6046; Sigma-Aldrich) and Nutrient Mixture F-12 
Ham (N4888; Sigma-Aldrich) supplemented with 10% fetal calf serum (Thermo Fisher Scientific), and 100 mg/mL 
penicillin/streptomycin, 37 °C, 5% CO . For experiments, cells were seeded on non-coated cover glasses and cultured 
in high-glucose Dulbecco Modified Eagle’s Medium (D1145; Sigma-Aldrich) supplemented with 10% fetal calf serum, 
2mM glutamine, and 100 mg/mL penicillin/streptomycin. Before fixation, cells were incubated for 2 h in the CDK1- 
inhibitor RO-3306 (10 uM in final concentration; SML0569; Sigma-Aldrich). 

Immunostaining. After 8 h, cells formed a closed monolayer. To increase the cell density within the monolayer and 
to attain a buckling instability, cells were cultured for a total of 25 h and then fixed for 15 min in 4% paraformaldehyde 
(43368; Alfa Aesar) in phosphate-buffered saline (PBS). After fixation, cells were permeabilized for 10 min with 0.1% 
Triton-X 100 and blocked for 60 min with 1% bovine serum albumin in PBS. E-cadherin was visualized using an E- 
cadherin rabbit antibody (1:500 ratio; 24E10; Cell Signalling) followed by staining with Alexa 532 antirabbit secondary 
antibody (1:500 ratio; A-11009; Invitrogen). F-actin was stained with Alexa Fluor 647-labeled phalloidin (1:500 ratio; 
A22287; Invitrogen) and the DNA with DAPI (1:1000 ratio; Sigma- Aldrich). 

Microscopy. Samples were imaged at high resolution on a home-build optical microscope setup based on an 
inverted Axiovert200 microscope body (Carl Zeiss, Oberkochen, Germany), a spinning disk unit (CSU-X1; Yokogawa 
Electric, Musashino, Tokyo, Japan), and an emCCD camera (iXon 897; Andor Labs, Morrisville, NC). 1Q-software 
(Andor Labs) was used for setup-control and data acquisition. Illumination was performed using fiber-coupling of 
different lasers (405 nm (CrystaLaser, Reno, NV), 514 nm (Cobolt AB, Solna, Sweden), and 642 nm (Spectra-Physics 
Excelsior; Spectra-Physics, Stahnsdorf, Germany)). Cells adhered on a cover glass were inspected with an EC Plan- 
NEOFLUAR 40 1.3 Oil Immersion Objective (Carl Zeiss). 

Image analysis. Cell segmentations and the height-to-radius profile analysis were performed using written scripts 
in Matlab2018a. Cell boundaries were identified from a maximum intensity projection of the F-actin signal of a 
confocal z-stack of the top part of the dome. The height profile was determined by the averaged intensity of the 
F-actin signal of a radius-dependent annulus area per plane of the z-stack. Fiji software was used for the orthogonal 
view of the dome. 3D reconstructions were done by ImarisViewer9.7.0 and z-directions were corrected for spherical 
aberration and axial distortion [1]. 


II. HEIGHT FUNCTION OF MDCK DOMES 


MDCK GII epithelial cells were growing on a coverslip. We observe that, after they formed a closed monolayer 
on the substrate, non-planar features, so called domes, developed due to the superelastic properties of the cells [2]. 
We imaged domes of different heights and diameters (see Fig. 1). Depending on the focal plane between the inner 
and outer area of the monolayer, cells changed their cell-cell contact length and the number of nearest neighbors. As 
described in the main text, this leads to the presence of topological defects in regions of high curvature (Fig. 1A and 
B). Furthermore, we sometimes found the additional disordered and folded structures on top of domes, see for example 
Fig. 1C. We measured the height of the buckled cell layer above the flat substrate. The height-to-radius profile of 
the outer area of the monolayer shows with increasing radius an increase in height, see Fig. 1D. We assume that the 
buckling instability and the formation of additional structures on top of the dome is caused by the cells continued 
division and growth even after they form a closed monolayer. The growth and division would cause the buildup of 
stresses in the monolayer that are relieved by the buckling. This is in agreement with e.g. Ref. [3]. We speculate 
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Figure 1. Domes in MDCK layers. Pictures of three domes of different sizes in panels (A)-(C) for which we measured the 
height of the cell monolayer above the flat substrate. The height of the outer area of domes as a function of the radius (distance 
from the center of the dome along the flat reference plane) is shown in panel (D). The color of each plot corresponds to the 
label in the panels (A)-(C), cyan to the dome in Fig. 1D in the main text. Segmented cells of domes shown in panel (A) and 
(B) illustrate the presence of topological defects near the top of the dome. For the three smallest domes, namely the functions 
colored blue [shown in panel (B)], green [shown in panel (A)], and cyan (shown in Fig 1D in the main text), we see that with 
increasing height at the center both the radius of the dome and the width of the plateau of constant height near the center are 
increasing. Furthermore, the height is a monotonically decreasing function of the radius. Instead, for the largest dome [red and 
shown in panel (D)] additional chaotic and folded structures were observed on top of the dome such that the height is maximal 
away from the origin. Red: F-actin, green: E-cadherin, blue: nuclei. 


that, consistently with the mechanism outlined in the main text, a net positive topological charge could facilitate the 
out-of-plane deformation of the monolayer, hence the formation of domes, in combination with other system-specific 
mechanisms, such as the injection of fluid under the cell layer which results in a focal detachment [2, 4]. Furthermore, 
it is possible that a difference in ion-concentration surrounding the monolayer between apical and basal side can even 
support the dome formation in areas where topological defects appear, see also Ref. [4]. 


To support this speculation, we have we analyzed eleven domes and counted the total topological charge in the 
central region, where the Gaussian curvature is maximal and positive. The topological charge of a cell having c; sides 
is conventionally defined as q; = 6 — ci. Thus, pentagonal cells (i.e. c; = 5) have topological charge qi = 1, whereas 
heptagonal cells (ie. c; = 7) has charge q; = —1. The total topological charge is then computed by summing the 
individual charges of all the cells in the central region of a dome, hereafter referred to as “region of interest” (ROI): 
Q = Vicror% (see Fig. 2A). We found that all domes in our sample feature a positive total topological defect charge, 
with the mean charge of the eleven domes being Qmean = 9.3 + 4.9 (mean + standard deviation). 


Furthermore, as shown in Fig. 2B we find that the total topological charge of a dome is strongly correlated with the 
total number of cells per dome (linear correlation coefficient r = 0.85). Such a linear relation, is consistent with the 
hypothesis that the topological charge of the cellular monolayer effectively screened by its Gaussian curvature, i.e. 


Qx dAK. (1) 


dome 


To illustrate this point, one can approximate the dome as a spherical cap of radius R, so that K = 1/R? and 
dA = dQ R?, with dQ the infinitesimal solid angle. The right hand side of Eq. (1) equates then the solid angle AQ 
spanned by the ROI: i.e. Q ~ AQ (see Fig. 2A). This, in turn, is proportional to the number of cells it encloses, given 
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Figure 2. Topological charge of domes. (A) Sketch of topological defects on a dome. Heptagons are blue and pentagons 
are red. (B) The total topological defect charge per dome increases linearly with the number of cells constituting the dome. 
(C) The cells contributing a positive topological charge are mainly pentagons. The height of each bar in the histogram is the 
averaged probability of finding a cell with n nearest neighbors in the eleven analysed domes. The error bar represents their 
standard deviations. 


that Nees = Aror/Acen = R?AQ/Acen. Thus 


Q x (4) Neells . (2) 


Fig. 2C shows instead the probability distribution of the number of sides c; of the cells in the ROI obtained from a 
sample of eleven domes. We see that a large majority of the cells are either pentagonal or hexagonal with almost half 
of the observed population being pentagonal (probability of 0.42 + 0.09 (mean + standard deviation)). Since, defects 
larger than 6—fold cells are underrepresented this leads to all the domes observed having a positive total topological 
charge. 


III. FORCE BALANCE 


In general, the force balance on the surface is given by Vio’ = —S°* [5, 6],where øf denotes the surface stress and 
=ext external forces applied to the system. We use bold letters to refer to vectors in R and latin indices to refer to 
surface coordinates on M. For example, the stress tensor can be written as a 2 x 3 matrix. It is possible to decompose 
the stress into its tangential and normal component, ot = ote; + ofn. In the absence of external forces, the force 
balance can then be written as 


Vio’ = (Vio + Koi) ej + (Viol, = o Kiz) n=0. (3) 


Here, we used Vin = Ki ej and Vjej = —K;jn. From the free energy of the main text we find, following Ref. [7-10] 


a 
the tangential and normal stress tensors 


of, = Jij — KBH Kij + KpH? gij +KF [gij XV?X — XViV;x + VixV5Xx| ; (4) 
and 


oni = KBV:H + Kp [2H = 2K;)x + 2(Ki = 29,’ H) Vix + krKÍXV;X : (5) 


To derive the equations in the main text we use the well-known geometric identity Kgi; = 2H Kij — Kink ey and the 
commutator [V;, V7]x = —KV;x. Using these together with the other stress tensors presented in the main text, we 
find the hydrodynamic equation for the momentum density and the normal force balance by projecting the general 
force balance onto the tangential and normal direction, respectively. We will present a more detailed derivation in a 
future work. 


IV. DERIVATION OF HEIGHT EQUATION 


In this section we derive the height equation, Eq. (5) of the main text, from the general equations describing the 
dynamics of the active surface, Eqs. (2) and (3) of the main text. We work in the Monge gauge where a height field 
h(r) is defined as described in the main text. A +1 defect is at the center of a disk of radius R and we assume that 
|Vh| < 1. As shown in the main text, we can find the explicit solutions for the velocity and director fields given 
in Eqs. (4) from the Eqs. (2) of the main text and we are thus now concerned with rewriting Eq. (3) of the main 
text. In the small-gradient approximation we have K ~ 0, H ~ —V7h/2, and the metric is just the identity such that 
covariant derivatives are equal to the flat derivatives. It is then straightforward to see that 


ft = —7V?h(r) + Evih) ; (6) 


Furthermore, using x(r) = —logr/R, we find 


peste rO?h(r) — 0,h(r) l (7) 


r3 


Lastly, for the curvature tensor coupled to the hydrodynamic and active stress tensor, we find the following. We have, 
Kyo") = (Vi V jh) [Pa — n (Viv + Vi0')] , (8) 


and 
r) d e a 
Kujo’ = a(ViV jh) (50 -p'pi) ' (9) 


where 6;; is the flat polar metric. For the first term on the right-hand side of both Eq. (8) and Eq. (9) we have 
(V'VIh)6,; = V?h. Furthermore, in polar coordinates V;Vjh = ô;ðjh — Tk ðkh, with T}; the Christoffel symbol 
associated with 6;;. This is non-zero only if i = j = r, in which case V7h = Oh, or if i = j = ọ, then V2h = rô,h. 
Similarly, we have Vivi = Ov) + ËT? v". With v, = 0 and vy = v,(r) we then find that this is non-trivial only if 
i = r and j = ọ such that Vv? = "v?. Thus, since V Voh = 0, the second term in Eq. (8) vanishes. Hence, as 
mentioned in the main text, the velocity does not enter explicitly in the final equation for h. For the last term we 
have: 


sin? € 


ppi ViV jh(r) = cos? e8?h(r) + O,h(r) , (10) 


r 


where e = + arccos(—1/A)/2 as found above. Therefore, adding Eqs. (6)-(9) together, we find 


a > h r-1 A+1 KB 
= fE + fi + Ky (oT +0) = (P, + 5-7) nt a, 2h 4 ph 
0= fè t+ fa + j (0 +o ) at 5 IEN oan ee 2A or Dw ay 2 


which can be rewritten to yield the height equation, Eq. (5) in the main text. 


V. MAPPING TO PHYSICAL UNITS 


The model parameters in lattice units used for simulations are kg = 0.008, Ag = 0.02, xr = 0.02, 8 = 0.03, M = 
0.1,f =1,A = 1.1, and 7 = 5/3. By following previous studies [11, 12], an approximate relation between simulation 
and physical units (for an active gel of cytoskeletal extracts) can be obtained using L = 1 ym as the length-scale, 
7 = 10 ms as the time-scale, and F = 1000 nN as the force-scale. A mapping of some relevant quantities is reported 
in Table I. 


VI. MOVIES 


Movie 1: 3D dome animation. 3D rotation and reconstruction of the same MDCK GII monolayer as in Fig. 
1D of the main text. Double scaling of the z-direction was chosen to ease the visualization of the buckling instability. 
The color code is as follows. Red: F-actin, green: E-cadherin, blue: nuclei. 


Model parameters Simulation units|Physical units 
Shear viscosity, 7 5/3 1.67 KPas 
Elastic constant, KF 0.02 100 nN 
Flow-alignment parameter, | 1.1 dimensionless 
Diffusion constant, D = Ma |0.001—0.015 (0.0087 — 0.0128 ppm?s~! 
Activity, a 0 — 0.01 (0 — 100) KPa 
Table I. Mapping of some relevant quantities between simulation units and physical units. 


Movie 2: Profile and tessellation. Left: top-to-bottom confocal z-stack of the same MDCK GII monolayer as 
in Fig.1D of the main text. At a distance of 15.3 um to the coverslip, cells are segmented and correlated according to 
their number of nearest neighbors. Top-right: cross section of the dome. The white lines indicate the inner, middle 
and outer area of the monolayer. Bottom-right: height-to-radius profile of the buckled monolayer starting from the 
center of the dome (see also Fig. 1). The position (red circle) follows a sigmoidal fit of the outer area and moves 
according to the z-stack animation in the left panel. Red: F-actin, green: E-cadherin, blue: nuclei. 

Movie 3: Oscillation between buckled states. Oscillations between a cuspidal configuration with negative 
curvature at the center and a smooth configuration. From an initially flat state the interface buckles. Note that 
while the initial protrusion is growing the defect is slightly off-center, leading to the whole protrusion moving on the 
xy-plane. (The z-direction corresponds to the interface normal in the initial state.) After some time (t > 380000) 
the interface starts to periodically oscillate between these two configurations. The oscillation is accompanied by the 
growth and shrinking of additional thin protrusions at the top. The white vectors denote the polarization field, while 
the color code refers to the local magnitude of the flow according to the color bar on the right-hand side. 

Movie 4: Droplet nucleation. The large extensile active stresses lead to the rapid growth of protrusions that 
cannot be counteracted by the elasticity of the interface. This results in a periodic breaking of the interface and the 
consequent nucleation of a droplet after each rupture. There are two +1 defects in the polarization field on each 
of the droplets. These deform under the straining action of the active liquid crystal on the surface and eventually 
dissolve due to Ostwald ripening. The white vectors denote the polarization field, while the color code refers to the 
local magnitude of the flow according to the color bar on the right-hand side. 

Movie 5: Chaotic dynamics. Shown is an example of the chaotic dynamics found at very high activity. The 
complex dynamics in this regime are characterized by chaotic deformations of the membrane with the consequent pro- 
liferation of many protrusions at the membrane. The active stresses cannot overcome the elastic forces of the interface 
and the protrusions break off from the membrane and elongate under the straining effect of bending deformations 
in the polarization pattern. The white vectors denote the polarization field, while the color code refers to the local 
magnitude of the flow according to the color bar on the left-hand side. 
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